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ABSTRACT 

As the number of discovered extrasolar planets has been increasing, diversity of planetary systems 
requires studies of new formation scenarios. It is important to study satellite formation in circum- 
planetary disks, which is often viewed as analogous to formation of rocky planets in protoplanetary 
disks. We investigated satellite formation from satellitesimals around giant planets through iV-body 
simulations that include gravitational interactions with a circumplanctary gas disk. Our main aim is 
to reproduce the observable properties of the Galilean satellites around Jupiter through numerical sim- 
ulations, as previous iV-body simulations have not expla ined the origin of th e resonant configuration. 
We performed accretion si mulations based on the work of [Sasaki et al.l (|2010l ). in which an inner cavity 
is added to the model of iCanup fc Ward] (|2002t 120061 ) . We found that several satellites are formed 
and captured in mutual mean motion resonances outside the disk inner edge and are stable after rapid 
disk gas dissipation, which explains the characteristics of the Galilean satellites. In addition, owing 
to the existence of the disk edge, a radial compositional gradient of the Galilean satellites can also be 
reproduced. An additional objective of this study is to discuss orbital properties of formed satellites 
for a wide range of conditions by considering large uncertainties in model parameters. Through nu- 
merical experiments and semianalytical arguments, we determined that if the inner edge of a disk is 
introduced, a Galilean-like configuration in which several satellites are captured into a 2:1 resonance 
outside the disk inner cavity is almost universal. In fact, such a configuration is produced even for 
a massive disk > 10 4 g cm~ 2 and rapid type I migration. This result implies the inevitability of a 
Galilean satellite formation in addition to providing theoretical predictions for extrasolar satellites. 
That is, we can predict a substantial number of exomoon systems in the 2:1 mean motion resonance 
close to their host planets awaiting discovery. 

Subject headings: planetary systems: formation - Planets and satellites: formation - Planet-disk 
interactions 



1. INTRODUCTION 

Origins of satellites around outer giant planets such 
as Jupiter and Saturn are important for determining the 
history of these objects, and their existence holds clues to 
the origins of such planets. Satellites around Jupiter and 
Saturn, in fact, are suitable targets for formation theo- 
ries because their orbital parameters and co mpositions 
are well documented (e.g., iSeidelmannl [19 92) compared 
to those of exoplanets. Therefore, we have developed a 
formation theory that compares physical parameters. 

Besides satellites in the solar system, those in extra- 
solar systems, known as exomoons, are an additional 
focus of this study. Exomoons have recently attracted 
great interest because of their habitability. The num- 
ber of detected extrasolar giant planets is much larger 
than rocky exoplanets, with the exception of Kepler can- 
didates, and some orbit within the habitable zone (HZ), 
i.e., an orbital region in which the stellar flux is suffi- 
cient to maintain liquid w ater on the surface of a planet 
(e.g., IKasting et al.l [19931) . If rocky satellites orbit gi- 
ant planets in the HZ, they are potentially habitable. 
According to studies on the possibility of life-bearing 
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moons (Wil liams etall 119971 lKalteneggerll2000D . satel- 
lites around extrasolar giant planets in the HZ might be 
habitable if their size is sufficient (> 0.2 M ffi ). An ob- 
vious exception is Titan, which has a dense atmosphere 
even though its mass is ~ 0.02 M®. 

A detaile d study on the detect ability of habitable ex- 
omoons by Kipping et al.l ([2009D determined that exo- 
moons in the HZ down to 0.2 M ffi may be detected by 
transit timing effects, including transit timing variations 
and transit duration variations, with the expected perfor- 
mance of the Kepler space telescope. In fact, an ongoing 
observational project uses photometry da ta of Kepler to 
detec t and investigate exomoon signals (Kipping et al. 
2012). In addition, the possibility of screening the at- 
mosphere of exomoons for habitabilit y with transmissio n 
spectroscopy has been discussed bv iKalteneggerl (|2010f l. 
who discovered that the number of transits needed to 
detect biomarkers on an Earth-like exomoon under ide- 
alized conditions and viewing geometry is feasible using 
the James Webb Space Telescope (JWST) for the sample 
of the closest M stars. Thus, it is important to predict 
the characteristics of satellites, such as mass, composi- 
tion, and orbital location, that are likely to exist around 
extrasolar giant planets. 



2 



Ogihara et al. 



Table 1 

Properties of the Galilean Satellites 



Satellites 


a (R P ) 


M (icr b M P ) 


P (g cm ") 


e 


i (rad) 


Io 


5.9 


4.7 


3.53 


0.0041 


0.0007 ! 


Europa 


9.4 


2.5 


2.99 


0.01 


0.0082 , 


Ganymede 


15.0 


7.8 


1.94 


0.0015 


0.0034 


Callisto 


26.4 


5.7 


1.83 


0.007 


0.0049 



Note. — Values are taken from [Schubert ct al. ( 2004) and I Yoderl 
(1995). a, M, p, e, and i are the semimajor axis, mass, material density, 
eccentricity, and inclination of the satellites, respectively. 



Regular satellites, which follow prograde and rela- 
tively close orbits up to tens of planetary radius with 
little orbital eccentricities and inclinations, are be- 
lieved to have formed as by-products of planets' for- 
mation within a circump lanetary accretion disk (e.g., 
ILunine fe Stevenson! 119821 ) . Circumplanetary disk mod- 
els currently fall into two categori es: (i) solids-enhanced 
minimum mass (SEMM ) model (Mo saueira fc Estradal 
200310; lEstrada et all I2009I ). in which a station- 
ary disk is assumed; and ( ii) gas-starved disk model 
HCanup fc Wardl[200a . [20(51 IWard fc Canud 120101) . also 
known as a slow inflow disk, in which an actively supplied 
gaseous disk is considered. 

In the SEMM model, a disk is composed of two dif- 
ferent zones: an optically thick inner region with a peak 
surface density ~ 10 5 g cm -2 and an optically thin outer 
region extended to a fraction of the Hill radius of the 
planet (~ 0.2 rjj). The transition from the inner disk 
to the outer disk is assumed to occur at the location be- 
tween Ganymede and Callisto for the Jovian disk. In this 
model, the disk is static, that is there is no inflow from 
circumstcllar orbits. The inner high density region with 
very low turbulent viscosity (a ~ 10~ 6 — 10~ 5 ) is prac- 
tically inviscid. Solid components of the disk are sup- 
plied by ablation and capture of planetesimal fragments 
passing through the disk. In the gas-starved disk model, 
however, a low-mass viscously evolving disk with a peak 
surface density of ~ 100 g cm -2 is considered. This disk 
is continuously supplied by the ongoing inflow of gas and 
dust particles from the protoplanetary disk. Disk condi- 
tions evolve in a quasi-steady state mode in response to 
the rate of gas accretion and turbulent diffusion, leading 
to lower disk gas densities than in the SEMM model. In 
this study, we adopt a fiducial model for gas surface den- 
sity based on the gas starved disk model (~ 100 g cm -2 ). 
Note that if low disk viscosity or rapid gas inflow is as- 
sumed, the gas surface density can become higher than 
100 g cm -2 (see Equations fl3J) and (gj). 

Our main focus in this paper is the formation of 
Galilean-like satellites, although we extend our theory 
to exomoons and also comment on the Saturnian system 
later in this paper. Physical properties of the Galilean 
satellites are shown in Table [TJ Orbital and composi- 
tional characteristics can provide valuable constraints on 
formation scenarios. The orbital resonances among the 
satellites Io, Europa, and Ganymede present a fascinat- 
ing dynamic system. Io-Europa and Europa-Ganymede 
are in the 2:1 mean motion resonance, which causes their 
successive conjunctions to occur near the same longitude. 
These satellites are also in the Laplace resonance, which 
denotes the 1:1 commensurability between the rates of 
motion of the Io-Europa and Europa-Ganymede conjunc- 



tions; thus, a triple conjunction never occurs. These res- 
onances suggest that the Galilean satellites underwent 
orbital migration either during or after their formation. 
With regard to composition, a progressive increase in 
"satellite ice mass fraction with increasing distance from 
Jupiter is observed. Io is composed of rock, and Europa, 
Ganymede, and Callisto are believed to contain approx- 
imately 8%, 45 %, and 56% ice and water by mass, re- 
spectively (e.g.. iSchubert et aLl l2004). In addition, Cal- 
listo's interior structure allows it to constrain its accre- 
tion timescale. Data from the Galileo spacecraft suggest 
that Callisto appears to be largely undifferentiated, pro- 
vided that the satellit e is in hydrostatic equilibrium (e.g., 
lAnderson et al.l I2001D . Therefore, Callisto must have 
avoided melting in its entire history. Estimates of the 
temperature increase associated with accretional heating 
show that Callisto could have remained unmelted during 
formation if it s accretion was on a t imescale longer than 
5 x 10 years (jBarr fc Canud 120081 ). In addition, eccen- 
tricities and inclinations are considerably small, which 
are considered to be damped by the surrounding gas neb- 
ula and/or tidal dissipation within the satellites. 

Formation of the Galilean sat ellites has been investi - 
gated from several perspectives. iCanup fc Ward! (2006) 
performed TV-body simulations of satellite accretion from 
satellitesimals in a starved disk that included the ef- 
fects of a gas disk. Because the type I decay and accre- 
tion timescales are shorter than the disk lifetime, several 
generations of satellite formation and migration are re- 
peated before gas disk dissipation. As the infall of gas 
and solid materials wanes because of global depletion of 
the circumstellar disk, the circumplanetary disk is de- 
pleted through viscous spreading and/or photoevapora- 
tion. Thus, the c urrent satellites are the survivors of the 
final generation. ICanup fc Ward! (|2006l ) also found that 
the total satellite mass fraction to the host planet is reg- 
ulated to ~ 10 -4 Mp (Mp; mass of the planet), which 
is only weakly dependent on model parameters and is 
consistent with Jovian and Saturnian systems. 

Several aut hors attempt to cl arify the origin of res - 
onant states. lYoderl (| 1979( 1 and lYoder fc Pealel (|1981( 1 
proposed that Io captured Europa into the 2:1 reso- 
nance through differential orbital expansion due to grav- 
itational tides raised on Jupiter. With continuing or- 
bital expansion, Europa eventually encountered the 2:1 
mean motion commensurability with Ganymede, result- 
ing in capture in to the Laplace resonance . In contrast, 
iGreenberd (|1987D and iPeale fc Leel (|2002( 1 advocated a 
primordial origin of the resonant states through differen- 
tial inward type I migration in the gaseous disk. Since 
the type I migration timescale is inversely proportional to 
the satellite mass, the largest satellite, Ganymede, under- 
went the most rapid drift and caught up with the inner 
satel lites, leading to c aptur e into the resonances. How- 
ever, ICanup fc War d (2006) performed detailed TV-body 
simulations, which is the first TV-body study of satellite 
accretion, and it is suggested that resonant states can- 
not be produced in their calc ulation^]. It seems that the 
initial orbital configuration of lPeale fc Led (|2002l ) cannot 



1 ICanu p & Ward ( 2006) did not analyze whether the satellites 
evolved in resonances or not, but we performed jV-b ody simulations 
under the same condition as Canup & Ward ( 2006) and confirmed 
that resonant relationships were hardly established. 
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be established during Canup & Ward's iV-body simula- 
tions, because when a Ganymede-mass satellite forms in 
outer region, satellites that reside in inner region do not 
grow to the mass comparable to those of Io and Europa. 

As shown in recen t studies on format io n of 
close-in super-Earths dTerauem fc Papaloizoul [20071 : 
lOgihara fc Idall2009t iPapaloizou fc Terqueml 12010 ) . mi- 
grating bodies can be tra pped in mean motion resonances 
near the disk inner edge. iSasaki et al.l (|2010l ) introduced 
the concept of the disk inner cavity in their investiga- 
tion of satellite formation and found that four or five 
satellites are generally formed by being trapped in reso- 
nances, whic h agre es well with the Galilean satellites. 
ISasaki et al.l (|2010D adopted a semianalytical method 
rather than direct orbital integration and performed 
Monte Carlo simulations. Although this approach pro- 
vides a powerful means for overall understanding of satel- 
lite formation, it is not adequate for effectively discussing 
detailed features of resonant trapping, because partic- 
ularly in the formation of Galilean satellites, similar- 
sized satellites interact with each another, and multi- 
ple bodies are eventually trapped in resonances. Several 
studies outlined the prediction of capture probability for 
resonances as a function of initial particle eccentricity 
in the adiabatic limit, where the migration timescale is 
much longer than the res onant libra tion tim e scale (e.g., 
Murray fc Dermott]ll999D . Recently. lOuillenl (|2006l ) and 



Mustill fc Wvattl (|2011l ) studied a rapid migration case 



in the restricted three-body problem in the limit of low 
initial particle eccentricity and semianalytically obtained 
capture probability using a Hamiltonian model. How- 
ever, it is uncertain whether their formulae can be used 
in a case in which bodies have comparable masses and 
the effect of eccentricity damping is included. There- 
fore, the TV-body simulation is the only currently avail- 
able method to accurately examine the detailed features 
of resonant trapping in systems in which multiple bodies 
with comparable masses interact with each other. 

Thus, we conduct TV-body s imulations of satellite ac- 
cretion in a gas-starved disk (jCanup fc Ward! 120021 ) by 
adding an inner cavity to the original model, to clarify 
the formation scenario with the disk edge. In particular, 
we assess the probability of primordial formation of the 
resonance relationship. Moreover, we attempt to repro- 
duce additional properties of Galilean satellites such as 
the total number and the compositional gradient. 

Furthermore, we examine the dependence of model pa- 
rameters on the properties of formed satellites. Recent 
studies on disk-planet interactions have suggested that 
torque exerted on a planet from a disk can be signif- 
icant ly altered because of disk str ucture dMasset et al 
2006) and radiative effects (e.g., iPaardekooper et al 
20101 ). resulting in a discrepancy with the type I migra- 



tion s peed predicted by linear calculations (jTanaka et al.1 
2002) . In an actively supplied disk, in which equilibrium 
between gas inflow and viscous diffusion is assumed, gas 
surface density depends on the inflow flux and gas vis- 
cosity. Given such uncertainties, we discuss satellite for- 
mation under a wide range of parameters. The results 
of this study are instrumental for effective prediction of 
orbital properties of extrasolar satellite systems. 

In Section [3] of this paper, we detail our model and 
numerical method; in Section[3j we present the results of 
iV-body simulations; in Section [4j we discuss the depen- 



dence of model parameters on the properties of formed 
satellites; in Section we briefly discuss satellite for- 
mation for the Saturnian system; and in Section [6l we 
summarize our conclusions and the implications of this 
study. 

2. MODEL AND CALCULATION METHOD 



ISasaki et al.l (|2010f ) considered two different models 
that correspond to Jovian and Saturnian systems: Sur- 
viving satellites around Jupiter were formed in a rela- 
tively massive protosatellite disk with an inner disk cav- 
ity, and those around Saturn were formed in a less mas- 
sive and cooler disk without an inner disk cavity. These 
disk conditions are inferred from gap openings along the 
host planet's orbit in the circumstellar protoplanetary 
disk and through observations of classical T-Tauri stars 
(CTTSs) and weak-line T-Tauri stars (WTTSs). 

Gap openings require that both viscous and t hermal 
conditions be satisfied ()Lin fc Papaloizou 11985') . The 
critical masses are expressed as ( Ida fc Lm1l2004l ): 



M, 



g,vm 



M, 



g,th 



where M g y ls and M 9 . t h are the viscous and the ther- 
mal conditions, respectively. Although these masses have 
some uncertainties, M g y ls and M g _ t h are larger at Sat- 
urn's orbits than at Jupiter's orbit. Since Saturn's mass 
is three times less than Jupiter's one, it is plausible that 
Jupiter opened a gap to halt its growth while Saturn did 
not and Saturn's growth wa s terminated by g lobal dis- 
sipation of the solar nebula (Sa saki et al.ll2010l ). As will 
be discussed in Section |2~T1 gap openings lead to rapid 
depletion of protosatellite disks; satellites formed in a rel- 
atively massive protosatellite disk could survive around 
Jupiter, while those formed in a late-stage less massive 
d isk could survive around Saturn. 

iHerbst fc Mundtl (|2005l ) observed that spin periods of 
young stars show bimodal distribution in peaks at about 
1 week and 1 day. They suggested that the stellar mag- 
netic field of the 1-week-period stars is coupled with the 
circumstellar disk and is sufficiently strong to transfer 
spin angular momentum to the disk and open a disk in- 
ner cavity, whereas the disks around the 1-day-periq d 
stars do not have the cavity (see also iHartmannl 12002') . 
They also suggested that the 1-week and 1-day stars may 
correspond to CTTSs and WTTSs, respectively, enabling 
massive CTTS disks with a cavity to evolve to less mas- 
sive WTTS disks without a cavity. In this evolution 
analogy, protosatellite disks may have a cavity in rela- 
tively early stages that disappears as the disk evolves, 
although this scenario includes a large uncertainty. Be- 
cause the current spin rate of Jupiter is much slower 
than its break-up spin rate, magnetic coupling between 
Jupiter and t he circum- Jovian disk is an obvious assump- 
tion. In fact. lTakata fc Stevenson! (|1996f ) showed that the 
circum- Jovian disk may have contained an inner cavity. 
The satellites formed in such a disk survived in case of 
Jupiter. The circum- Saturnian disk may once have con- 
tained a similar disk in early stages. However, if the disk 
evolved to a less massive disk without a cavity, the satel- 
lites formed in the massive disk would have fallen onto 
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Saturn, and the surviving satellites would have formed 
in the late-stage disk. We comment on the formation of 
the Saturnian system in Section In the rest of this 
paper, we focus on the model for Jupiter. 

2.1. Disk Model 
2.1.1. Gas Surface Density and Disk Temperature 

As stated in Section [TJ we adopted the actively sup- 
plied gas-starved disk model in o ur calculation s . Ou r 
disk model is based on that by Sasa ki et al.1 (|2010f ). 
in which the disk i nner c avity is added to the model 
by iCanup &; Ward! (|2002| ) with other slight modifica- 
tio ns. In this section , we briefly summarize the model 
bv lSasakiet^(l20Trl . 

For the case in which gas inflow from the circum- 
ste llar disk is limi t ed to the region between r ln and 
r c , ICanup fc Ward! (|2009t ) suggested that the inner and 
outer radial boundaries of the inflow region for the proto- 
Jovian disk are 4 Rj and 34 Rj, respectively. This the- 
ory is consistent wit h the th r ee-dim ensional hydrody- 
namic simulation by iMachidal (|2009f ). i.e., r c ~ 22 Rj, 
wh ere Rj is the physical radius of Jupiter. As ad opted 
by ICanup fc Wardl (pOOfih and ISasaki et all ffiwh . we 
used r c = 30 Rp in this study. The total inflow rate 
is expressed as Fp = Mp/tq, where tq is the gas 
inflow timescale; then the infall flux per unit area is 
Fn = Fp/ir{r1 — rf n ) ~ Fp/nr^. Here we ignore the 
radial dependence of the inflow because the resolution of 
the current hydrodynamic simulations is insufficient to 
constrain the radial dependence of the infall, as will be 
discussed in Section [5] 

A disk is formed with the inflowing gas, which diffuses 
viscously. If the viscous diffusion timescale is shorter 
than the characteristic timescale over which the inflow 
changes, the gas disk can be described as a steady accre- 
tion disk, in which the inflow and viscous diffusion are 
equilibrate d. According to the st eady-state disk model 
derived by ICanup &; Wardl ((20021 ) . the gas surface den- 
sity of the disk is approximately given by (Appendix IA"|) 
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where v, M, r, and R are the viscosity, mass, radial dis- 
tance from the planet, and physical radius of the bodies, 
respectively. The subscripts "P" and "J" denote quanti- 
ties of the host planet and Jupiter, respectively. Because 
the magnitudes of turbulent viscosity and the gas inflow 
rate are not well determined, we introduce the scaling 
factor / g for gas surface de nsity. We adopt an alpha 
model for the disk viscosity (Shakura & Sunyacv 11973T) 
v = ac s H ~ ac 2 s /^k, where H and fix are the disk 
scale height and the Keplerian angular velocity, respec- 
tively, and the sound velocity c s is derived thorough the 
temperature distribution T. T is determined by the bal- 
ance between viscous heating and blackbody radiation 



(Appendix [A"| as follows: 
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2.1.2. Gas Dissipation and Disk Inner Edge 



Gas surface density depends on the inflow rate; thus, 
the decay of the inflow leads to a decrease in gas sur- 
face density (and disk temperature). In our calculations, 
we generally assume an exponential decay of the inflow 
rate as Fp x exp(— i/rd cp ), where the decay timescale 
Tdcp is set to be 10 6 years. We therefore adopt an ex- 
ponential decay of gas surface density with a dissipation 
timescale of 10 6 years. Altough more accurate depen- 
dence of the inflow rate on the surface density is repre- 
sented by Sg oc Fpjv cc F^ A , we used the same decay 
timescale for simplicity. 

Although infall flux decay associated with global dis- 
sipation of a protoplanetary disk is ~ 10 6 years, we in- 
troduce a significantly more rapid decay associated with 
inflow truncation by the gap opening in the protoplan- 
etary disk. After the gap opening, the protosatcllitc 
disk would be rapidly depleted on the viscous diffusion 
timescale r^is, 



Tdiff 



10' 



Go- 3 ) 



yr- 



(6) 



The gap opening timesc ale itself would be comparable to 
Tdiff ((Sasaki et al.ll2010l) ; therefore, both timescales are 
much shorter than accretion and migration timescales of 
satellites. In our calculations, the disk gas is abruptly 
dissipated with Tdiff, such that the satellites are "frozen" 
at that time. Note that although hydrodynamic simula- 
tion suggests that the gap ope ning could not comp letely 
truncate the infalling gas (e.g.. iLubow et al]|1999ft . even 
in this situation, results that are shown in this study 
remain unchanged because the infall rate would be re- 
duced by several orders of magnitude. As long as the 
inflow cutoff timescale is shorter than accretion and mi- 
gration timescales, which are ~ 10 5 years, the results are 
not affected. That is, the results are not sensitive to the 
decay rate of the inflow. 

As the location of the inner edge of the protosatcllitc 
disk is not well constrained, in our calculations, we set 
the disk inner edge at 5 Rp, which is slightly outside 
the current corotation radius of Jupiter (~ 2.25 Rj). 
At the disk inner edge, the gas surface density of the 
disk smoothly vanishes with a hyperbolic tangent func- 
tion with the width Ar. We adopt Ar = 0.2 Rp, which 
is comparable to the disk scale height H. In addition, 
satellites that migrate inside 4 Rp are disregarded from 
calculation mainly for computational reasons, which is 
also described in Section [2. 4. 21 

2.2. Solid Inflow Model 

Based on the model used by ICanup fc Wardl ((2006) i 
solid bodies are added to the calculation between r; n = 
5 Rp and r c = 30 Rp at a rate -Fpryico/d. in/100, where 
/d,i n and 77i CG are the scaling factors of the amount of 
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Table 2 

Parameters for Each Simulation 



Model (runs) 


k 


/d,in 


C'i 


ice Line 


model l(a,b,c) 


1 


1 


1 


No 


model 2(a,b,c) 


1 


1 


1 


Yes 


model 3(a,b,c) 


10 


1 


1 


No 


model 4(a,b,c) 


1 


0.5 


1 


No 


model 5(a,b,c) 


1 


2 


1 


No 


model 6(a,b,c) 


1 


1 


0.1 


No 


model 7(a,b,c) 


1 


1 


10 


No 



Note. — Three runs referred to as a, b, 
and c (thus, a total of 21 runs) are performed 
for each model. / g and /d,i n are scaling fac- 
tors for gas surface density and the solid in- 
flow, respectively. C\ is the scaling factor 
for the type I migration speed. For model 2, 
solid enhancement outside the ice line is con- 
sidered. 



solid material and the increase of solid materials due to 
ice condensation outside the "ice line," respectively. The 
solid infall flux per unit surface area is assumed to be 
uniform because we assume a uniform infall flux of gas 
Fp, which is stated in Section [2.1.11 When ?7ice/d,m = 
1, the gas-to-solid ratio in the inflow is assumed to be 
100. However, the "effective" gas-to-solid ratio in the 
disk is not fixed to 100, because gas surface density is 
also parameterized by / g and the solid materials which 
are decoupled from the gas evolve independently, (see 
Tablc[2]for the amount of gas of each run.) We assume an 
exponential decay of the solid inflow with the timescale 

Of Tdep- 

Although the infall of solid components is the contin- 
uous dust infall in the gas flow, it is mimicked by adding 
bodies with mass M ad d = 5 x 10 _7 ryi CO (r/15i?p) 3 M P 
with a random azimuthal location at s ome time inter- 
val ( see Supplementary Information in iCanup fe Wardl 
12006( 1. Note that although a uniform solid inflow flux 
per area is assumed, radially dependent masses are used 
for adding bodies to avoid an imbalance in the number 
of adding bodies between radial zones. A total of 2,000- 
6,000 bodies are calculated for single runs, depending on 
the runs. The initial velocity dispersion of the proto- 
satellites is set to be their escape velocity. The corre- 
sponding initial eccentricity e and inclination i are given 
by v csc = Ve 2 + i 2 v K with e = 2i (|Ida fc Makinolll992[) . 

In some runs, solid enhancement outside the ice line 
is adopted. Assuming the condensation temperature of 
~ 160 K, the location of the ice line is derived from Equa- 
tion ([5]), and we consider that the ice line moves inward 

as the gas inflow rate decreases. Because Td cx F^ 4 and 



Td cx r 3 / 4 , the location of the ice line is proportional to 
30 exp ( - 1 /3r d ep ) Rp ■ 



F i/3. 

-Tp . / 



2.3. Orbital Integration 

The orbits of satellitcsimals are calculated by numeri- 
cally integrating the equation of motion of the particle k 
at rk in planetocentric coordinates, 



d 2 rk rk 



ig + -^tidc 1 



rk 



(J) 



where k,j = 1, 2, the first term on the right-hand 
side is the gravitational force of the central planet, the 
second term is mutual gravity between the bodies, and 
the third is an indirect term. .Fdamp, -Fmig, and .Ftide are 
specific forces due to gravitational eccentricity damping, 
semimajor axis damping (type I migration), and tidal ec- 
centricity damping, respectively, which will be explained 
in Section l2~4l 

For numeri cal integration, we use the fourth-order Her- 
mite scheme (jMakino fc Aar scth 19921) with a hierarchi- 
cal individual time step (jMakinol 119911 ). When physical 
radii of two spherical bodies overlap, we consider the 
bodies to merge, conserving total mass and momentum 
and assuming perfect accretion. The physical radius of a 
body is determined by its mass M and internal density 
p as 



',4tt p 



1/3 



(8) 



where we adopt p — 3 g cm 3 . 

2.4. Characteristic Timescales 
2.4.1. Gravitational Interaction with Disk Gas 

We consider damping of orbital eccentricity, inclina- 
tion, and semimajor axis due to disk-satellite interac- 
tions. A satellite gravitationally perturbs the disk gas 
and excites de nsity waves, which damp e,i, and o of the 
satellite (e.g-IGoldreich fc TremaineTfT98a IWardl[l98l 
IArtvmowiczlll993l) . 

The force formulae for e-damping and i-damping 
(Frfamn) are the sam e as that in Equations (2)- (4) in 
Oalh ara et al.l (120101) . and the e-damping timescale is 
(jTanaka fc WardE 004') 
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The force formula for q- damping (F m j^ ) is th e same as 
that in Equation (3) i n lOgihara et all (120101) . and the 
migration timescale is (jTanaka et al.ll2002T ) 
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where S g Cv r 



5 x 10 6 yr 

~ q and q — 3/4 is adopted from Equa- 
tion (|3|). C\ is a scaling factor to allow for retardation 

2 This damping is sometimes called "tidal damping." Be- 
cause damping of orbital elements through tidal energy dissipation 
within planets and tidal interactions between planets and satellites 
through their tidal deformation are also considered in this study, 
we use the term "tidal" only when we refer to the tidal energy dis- 
sipation in planets and the tidal interaction between planets and 
satellites. 
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and acceleration of type I migration. We include neither 
reverse torque of t ype I migration near the disk edge 
(|Masset et al.ll2006l) nor the radiative effect through the 
entropy gradient (jPaardekooper et al.ll201Ct ); such issues 
will be resolved in future studies. Instead, we handle t a 
as a parameter in some runs and focus on the effect of the 
eccentricity trap, which is a physical mec hanism to halt 
type I migration near the disk inner edge (Ogihara ct al. 
I2010T ). This trapping arises when a body with an ellip- 
tical orbit straddles the disk inner edge; the body gains 
a net positive torque from the gas disk, which can com- 
pensate for negative type I migration torques that are 
exerted on a reso nantly interact i ng convoy of bodies. 

As d iscussed in lQgihara et all (|2010l ). lTanaka fc Ward! 
(|2004f ) assumed small e and i (<C Cs/vk) to obtai n 
the damping rate of eccentricity. lOstrikerl (I1999D . 
iPapaloizou fc Larwoodl <j2000f l. and iMuto et all (|2011[ ) 
showed that the formulae of e-damping and a-damping 
must be multiplied by correction factors in the case of 
supersonic flow. When the correction factors are added 
to t e and t a , it is suggested that the eccentricity trap 
becomes more efficient because the correction factor for 
t a is a stronger function of ev\^/c s than that for t e . How- 
ever, we neglect these factors for simplicity. 

2.4.2. Tidal Dissipation 

Tidal dissipation within satellites plays an important 
role in satellite evolution. Thus, we incorporate tidal 
eccent ricity damping of satellites a s a force acting on 
them (Papaloi zou fc Ter quem 20101): 
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wher e the damping timescale is (jGoldreich fc Soteil 
1966) 
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Here, jl (~ 2>j2ki; k 2 is the love number) and Q are the 
effective rigidity and the tidal dissipation function of the 
satellites, respectively. Although there is an uncertainty 
in Q, it is likely tha t both Q and jl are of the order of 
10-100 (|Yoderill995D . We thus use Q' = jlQ = 1000 for 
our JV-body simulations. 

Note that tidal torque on the satellite changes 
its semima.jor axis. The migration timescale is 
(|Goldreich fc Sotenfl966T) 
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where Qp and k 2 ^p are the tidal dissipation function and 
the love number of the central planet, respectively. The 
direction of migration is determined by the position of 
the satellite with respect to the corotation radius with 
the spin of the host planet. For example, if the satel- 
lite is on the inside of the corotation radius, the satellite 
moves inward. Because the migration timescale is gen- 
erally much longer than the simulation time, we neglect 
this effect. The change in semimajor axis due to dissi- 
pation within the satellite is also ignored because this 
timescale is also long except for very close-in satellites. 
The migration timescale has a strong dependence on the 
semimajor axis; thus, extremely close satellites inside the 
corotation radius would fall onto the planet in a short 
time. As previously described, the current corotation ra- 
dius of Jupiter is ~ 2.25 Rp; satellites that migrate inside 
of this radius would be quickly lost to the planet. Be- 
cause of both computational cost and uncertainty in the 
corotation radius, we disregard bodies that reach 4 Rp 
from our calculation. 

3. TYPICAL RESULTS OF iV-BODY SIMULATIONS 

In this section, we show typical results with our fidu- 
cial parameters; the parameters for each simulation are 
listed in Table [2j First, we perform iV-body simulations 
in a gas disk that dissipates with the global dissipation 
timescale of r^ep = 10 6 years, assuming that the solid in- 
flow flux also decreases with the same timescale. Then, 
the host planet opens a gap in the circumstellar disk at a 
certain time to enable the gas and solid inflow fluxes to 
rapidly decrease with a disk viscous diffusion timescale 
of Tdiff ~ 10 3 years. However, because the timing can- 
not be determined, we randomly choose the moment and 
perform the subsequent simulation. We mainly focus on 
the Af-body results in a gas disk because we found that 
orbital configurations do not change after rapid gas de- 
pletion. 

3.1. Growth and Migration before Gap Opening 

Orbital evolution of satellites embedded in a gas disk is 
presented in this subsection. In our fiducial runs (models 
la, lb, and lc), the gas scaling factor / g = 1 is used, the 
solid inflow rate is fixed at /d,i n = 1, and the migration 
rate C\ = 1 is assumed to be that predicted by the linear 
theory. We note that these parameters are not unique 
values but contain some uncertainties. In Section 21 we 
will discuss the dependency of the results on these pa- 
rameters. Other effects such as the supersonic correction 
factor of gravitational damping formulae and solid en- 
hancement outside the ice line are not considered. 

The orbital evolution for model la is shown in the bot- 
tom panel of Figure [T] In this figure, 30 most massive 
satellites at each time are plotted as circles, the radii 
of which are proportional to the radii of the bodies. 
Tk(— 0.03 years) is the orbital period at a — 20 Rj 
around a Jovian mass planet (M = Mj). The dashed 
line at 5 Rp represents the location of the disk inner 
edge. When a satellite with mass greater than 10 -5 Mp 
falls onto the host planet, it migrates inside 4 Rp; the 
descent is indicated by an arrow. The evolution of mass 
is shown in the top panel of Figure [1] which can be useful 
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Figure 1. Top: Evolution of the solid mass for model la. The 
dashed line shows the total mass added to the calculation and the 
solid line shows the mass remaining within the disk. Bottom: Or- 
bital evolution of the satellites for model la. For each time, 30 
most massive satellites are plotted. The circles represent bodies, 
and the radii of the circles are proportional to the physical radii of 
the bodies. In the electronic version, bodies with M > 10 _6 Mp, 
M > 10~ 5 M P , and M > 10 _4 M P are expressed with blue, red, 
and green circles, respectively. The arrow indicates that the satel- 
lite (> 10 -5 M P ) migrates inside 4 R P . 

for physical understanding of what we observe here. The 
dashed line refers to the total mass added, and the solid 
line to the solid mass remaining within the disk. 

Although many satellitesimals are not initially placed, 
solid materials are continuously added to the calcula- 
tion, and the growth of satellitesimals proceeds in a man- 
ner similar to that o f planetesimals around a star (e.g., 
iKokubo fc Mai 11998ft . Through orbital repulsion of oli- 
garchs that form from planetesimals, their orbital sepa- 
rations are kept wider than five mutual Hill radii. The 
typical orbital separation is Ar ~ 10 — 15 rn, where Th 
is the mutual H i ll rad ius of oligarchs. As is discussed in 
ICanup fcW ard (20(H), it is important to note that owing 
to the low solid inflow rate, the growth rate of a satel- 
lite is governed by the solid inflow flux. The accretion 
timescale is given by 



Tacc,inflow^l-8 X lO^J / d ^ 
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growth timescale are presented in Appendix [B] 

When the satellite grows to a critical mass, which 
is determined by a balance between the accretion time 
Tacc, inflow and migration time t a , it begins to undergo or- 
bital decay. The critical mass M cr it is estimated from 
Equations (H} and (JT9]) as 
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More detailed descriptions of satellite growth and the 



The critical mass in our fiducial case is ~ 8 x 10 -5 Mp, 
which is comparable to the masses of of the Galilean 
satellites. Because satellites grow to some extent during 
migration, the factor (fyicc/d.in/Ci/g) in Equation (|20[) . 
which places a lower limit on the satellite mass, would 
be smaller than unity for the case of the Jovian satel- 
lite formation. A discussion of parameter dependence is 
given in Section [4j 

The migra ting satellites exhibit dynamics different 
from those in ICanup &: Wardl (|2006l ) because of the exis- 
tence of the disk inner cavity. The satellite that reaches 
the disk inner edge no longer experiences gas drag, and 
migration ceases. The subsequently migrating satellite is 
captured into a 2:1 mean motion resonance with the in- 
ner satellite, leading to excitation of eccentricities (~ 0.1) 
of both bodies. The inner satellite that resides near the 
disk inner edge gains a positive torque from the disk, and 
eventually, the satellites captured into the 2:1 mean mo- 
tion resonance remain in the disk. This trapping mech- 
anism of bodies near the disk edge has been both nu- 
merica lly and analytically investigated by lOgihara et al.l 
(2010). They refer to the positive torque gained by the 
satellite and the trapping mechanism as the edge torque 
and the eccentricity trap, respectively. It should be em- 
phasized that we do not establish a rigid wall at the disk 
inner edge or manually impose an edge torque to the 
satellite at the edge. We apply only instantaneous drag 
forces t hat have been physica lly examined in previous 
studies ( jTanaka fc Wardl 120041 ) onto all bodies. If the 
body at the disk edge has a nonzero eccentricity, it moves 
back and forth between the inner cavity and the gas disk. 
When it moves into the gas disk near the apocenter, the 
tangential velocity of the body is slower than the local gas 
velocity; therefore, the body suffers a tailwind. On the 
other hand, when it enters into the cavity near the peri- 
center, where the tangential velocity of the body is faster 
than the local circular Keplerian velocity, the body does 
not feel gas dr ag, which resu l ts in a net positive torque 
on the body. lOgihara et al.l (|2010ft determined that if 
the edge torque balances with the negative type I torque 
exerted on the bodies, the eccentricity trap works effec- 
tively and other migrating bodies can even be trapped 
without pushing the innermost body into the cavity. Fig- 
ure [T] shows that the eccentricity trap is retained, except 
at t ~ 7 x 10 6 Ik- In fact, this trapping ordinarily also 
occurs in other simulations not shown here. However, 
an upper limit exists on the total number of satellites 
that can be retained in the disk, which is discussed in 
Section 

We also found that the 2:1 mean motion commensura- 
bility is universal in our TV-body simulations. Figure [2] 
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Figure 2. Evolution of the period ratios for model la. Circles 
and triangles express period ratios of the innermost and the second 
innermost pairs, respectively. 



shows the evolution of period ratios. Open circles and 
triangles represent period ratios of the innermost and 
the second innermost pairs, respectively. Although some 
satellites are pushed inside the 2:1 resonant location by 
outer migrating satellites and are captured into closer 
first-order resonances (e.g., 3:2 and 4:3), almost all pairs, 
whose period ratios decrease above two, are locked in the 
2:1 resonance. Hence, the period ratio is two. However, 
the resonant angles do not always librate about fixed 
values but sometimes exhibit circulation. The commen- 
surate value and the occurrence of the eccentricity trap, 
which regulates the number of trapped bodies outside 
the disk edge, depends on eccentricity damping and mi- 
gration rates, which is discussed in Section [4] 

Satellites with mass M cr ;t begin to move inward and 
continue to grow during migration even after they are 
captured into a resonance. At 3 x 10 7 Ik when we 
stopped calculations, the masses were 2.4 x 10~ 4 Mp 
(innermost), 2.5 x 10~ 4 Mp (second innermost), 2.9 x 
1CT 4 M P (third innermost), and 3.0 x 10" 4 M P (fourth 
innermost), which are heavier than the current Galilean 
satellites. For satellites with masses comparable to those 
of the Galilean satellites to be formed, the gap opening 
time, in which the solid inflow is terminated, should be 
earlier than 3 x 10 7 Tk- In addition, the critical mass for 
migration (Equation (|2Tj|) ) must not exceed the masses 
of the Galilean satellites. A discussion of satellite mass 
appears in Section 14.31 

We performed two additional runs (models lb and lc) 
under the same conditions but with different random 
numbers for additional satellitesimal locations, and we 
confirmed that the overall evolutionary tendency is iden- 
tical. The averaged orbital distributions of the four in- 
nermost satellites (> 10 -5 Mp) at each time are plotted 
in Figure |3l Circles, triangles, and squares represent av- 
eraged values over three simulations at 1 x 10 7 , 1.5 x 10 7 , 



and 2 x 10 7 Tk, respectively. Error bars indicate la dis- 
persion. The semimajor axis distributions that are ob- 
tained at given times are similar to those of the Galilean 
satellites. The three inner satellites are usually in the 
2:1 resonance, and the fourth one is often captured into 
the 2:1 resonance. We set the disk inner edge at 5 Rp 
so that the innermost one is located at 5 Rp. This orbit 
is slightly interior compared with that of Io; thus, satel- 
lites in Figure [3] generally exist closer than the Galilean 
moons. The inclinations are kept small compared to the 
excited eccentricities due to resonant interactions, and 
they are comparable to the current inclinations of the 
Galilean satellites. Table [T] shows e and i of the Galilean 
satellites. 
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Figure 3. Averaged semimajor axes, eccentricities, and inclina- 
tions of the four innermost satellites at t = 1 X 10 7 Tk (blue circles) , 
t = 1.5 X 10 7 Tk (red triangles), and t = 2 X 10 7 Tk (black squares) 
for models la, lb, and lc, respectively. Error bars denote the 1<t 
dispersion. (A color version of this figure is available in the online 
journal.) 



3.2. Orbital Evolution after Gap Opening 

The inflow is abruptly cut off because of the gap open- 
ing; thus, the disk gas is depleted at a certain time dur- 
ing the evolution in Figure [1] Because the time of the 
gap opening is unknown, we randomly choose the mo- 
ment and continue calculations with gas dissipation on 
the timescale of ram = 10 3 years ~3x 10 4 Tk- 

Figure U shows subsequent evolution of model la in 
which the orbital calculation is restarted at 1.3 x 10 7 Tk- 
Considering the formation of the Galilean moons in 
which the fourth satellite, Callisto, is not in a mean 
motion resonance, it is likely that this satellite under- 
went little or no migration. To achieve little signif- 
icant migration, the gap opening timing, after which 
the circumplanetary disk dissipates on the timescale of 
Tdiff, must be comparable to the timing of the comple- 
tion of Callisto formation. We assume the time to be 



SATELLITE FORMATION AROUND GIANT PLANETS 



9 




5xl0 7 10 8 1.5xl0 8 5xl0 7 10 8 1.5x10' 

Time (T K ) Time (T K ) 

Figure 4. Subsequent evolution of model la with rapid gas dis- 
sipation. Left: Evolution of semimajor axis. Right: Evolution of 
eccentricities of the three inner satellites. The eccentricities of the 
innermost, the second innermost, and the third innermost satellites 
are expressed with the solid red line, the dotted blue line, and the 
dashed black line, respectively. (A color version of this figure is 
available in the online journal.) 

1.3 x 10 7 Tk- An alternative scenario is that Callisto 
was accreted from the slowl y inflowing materia ls dur- 
ing an imperfect gap opening (jSasaki et al.H20ldh . which 
naturally e xplains the inferred undifferentiated interior 
of Callisto (jBarr fc Canudl2008h . Both scenarios explain 
the non-resonant Callisto. However, because the timing 
of gap opening cannot theoretically be predicted, more 
detailed discussion of characteristics of Callisto is left for 
future study. 

The left panel of Figure |4] shows the semimajor axis 
evolution. Because the satellites are in mean motion res- 
onances and their orbital separations are wide (~ 15 rjj), 
orbital configurations hardly change. Although the semi- 
major axis can be slightly decreased because of tidal dis- 
sipation, such an effect is negligible. Tidal torque can 
also change the semimajor axis; however, we ignore this 
effect, as previously stated. The right panel shows the ec- 
centricity evolutions of the innermost (solid line), the sec- 
ond innermost (dotted line), and third innermost (dashed 
line) satellites. The e-damping timescale of the inner- 
most satellite, which is located at 5 i?p, is estimated as 
~ 10 6 Tk (Equation ([To])). This result is not strictly con- 
sistent with the actual damping time (~ 10 7 Tk) of the 
innermost satellite observed in Figure [4j because eccen- 
tricities of outer satellites in resonances are also dragged 
down by tidal dissipation in the innermost satellite via 
resonant interactions between them. Thus, the eccen- 
tricities of the resonant satellites become < 0.01, which 
is consistent with thos of the Galilean satellites (Table[T]). 
Regarding the evolution of the resonant angles, as the or- 
bital energy is dissipated because of tidal e-damping, the 
satellites evolve deeper into the resonances, leading to a 
small amplitude libration around 0° or 180°. Again, not 
all of resonant angles librate about fixed values. In three 
runs of all fifteen simulations (models 1, 3-6), the Laplace 
relationship, which indicates that 6*5 librates, can be ob- 



served. The Galilean-like configuration (61, 9 2 , O3, and 
#5 librate while 64 circulates) is reproduced by one run 
of all. The definition of the resonant angle is presented 
in Appendix [Cj 

3.3. Origin of Compositional Gradient 

By tracking the orbital evolution of all the satellites 
in the ./V-body simulation, satellite composition can be 
discussed. As mentioned in Section [TJ a compositional 
gradient exists in the Galilean satellites. That is, bulk 
densities of Io, Europa, Ganymede, and Callisto are 3.53, 
2.99, 1.94, and 1.83 g cm" 3 , respectively. The composi- 
tion of satellites is considered to be a reflection of disk 
temperature; thus, we discuss in this subsection the ori- 
gin of the compositional gradient in terms of disk temper- 
ature. Disks can be classified into three classes: (i) hot, 
in which only rocky materials exist; (ii) cool, in which the 
ice condensation line is located in the satellite-forming 
region, and rocky and icy objects coexist; (iii) cold, in 
which water vapor can condense into ice in all parts of 
the disk. 
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Figure 5. Evolution of semimajor axis and water mass fraction 
for model 2a. The color of each object corresponds to its water 
mass fraction. The long dashed line represents the location of the 
ice line rj ce . 

To observe the effect of the ice line, we performed 
an additional set of iV-body simulations (model 2) by 
considering the enhancement of solid materials outside 
the ice line assuming 77i CC = 3, and we track the accu- 
mulation and migration of satellites. As noted in Sec- 
tion I2.2[ we assume that the ice line moves inward ac- 
cording to r; cc = 30 exp(— i/3rd ep )-Rp, which indicates 
that a cool disk is considered. Figure [S] shows the evolu- 
tion of model 2a, in which the colors corresponds to the 
water (ice) mass fraction of satellites, from red (dry) to 
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light blue (up to 67%). (The colors are visible in the on- 
line version of the paper.) The long dashed line denotes 
the location of the ice line. The water content is calcu- 
lated using the following prescription for components of 
satellites that originated at r: 

Mwatcr _ ??icc - 1 Jo (r < r icc ) , , 

M ~ ??icc - \0.67(r >r icc ). V L > 

In the case without the disk inner edge, icy materials can 
reach the innermost region with the timescale of type I 
migration (~ 10 6 7k), and fine-tuning of timing is re- 
quired for the final survival of both rocky and icy satel- 
lites. The type I decay time should be comparable to 
the migration time of the ice line in the final state. How- 
ever, as shown in Figure [SJ icy satellites stay outer region 
it = 1 — 1.5 x 10 7 7k). Because satellites are lined up 
outside the disk inner edge, migrating icy satellitesimals 
are consumed by the outer satellites captured into mean 
motion resonances, leading to the deficiency of icy mate- 
rials in the inner satellites. Thus, contamination of inner 
satellites by icy materials is prevented owing to the disk 
inner edge, and we do not have to rely on the fine-tuning 
of the timing. We find a monotonic increase in the water 
content with an increase in the distance from the planet. 
Between 8 x 10 6 and 1.5 x 10 7 7k, the water mass frac- 
tion of the innermost satellite is ~ 5%, while that of the 
outer satellite is ~ 60%. This configuration is approx- 
imately consistent with the those of Galilean satellites. 
However, after 1.8 x 10 7 7k, the rocky inner satellites are 
lost to the planet because of the violation of the eccen- 
tricity trap. As a result, the surviving satellites contain 
a large quantity of ice (> 40%). 

Rather than the effect of the ice condensation line, we 
find in a companion paper (Nimmo, Ogihara & Ida, in 
preparation) that the impact erosion of ices during ac- 
cretion can produce the compositional gradient even in a 
cold disk. We find that the impact velocity increases with 
a decrease in the radial distance and exceeds 10 km s _1 
at the orbits of Io and Europa. This observation is at- 
tributed to the high Keplerian orbital velocity (cx r -1 / 2 ) 
and a large velocity dispersion. Thus, nearly all ice in the 
innermost body, which corresponds to Io, and most of the 
ice in the second innermost body, which corresponds to 
Europa, evaporates through collisions. The amount of 
ice vaporization depends on the mass ratio of the target 
and th e impactor, the impact velocity, and the impac t 
angle (|Kraus et all 120111: iNimmo fc Korvcanskvl 1201 21 ). 
Detailed analysis is presented in the companion paper, 
in which the evaporation of ice owing to tidal heating is 
also discussed. 

To summarize, to reproduce the compositional gradi- 
ent by considering ice condensation outside the ice line, 
which indicates a cool disk, the following two conditions 
are necessary. First, the location of the ice line, in which 
there is an uncertainty, should be sufficiently far to en- 
able the formation of a large amount of rocky satellites- 
imals. Second, once icy materials move into the inner 
satellite formation region, the eccentricity trap of a reso- 
nant convoy should not be broken. We also find that fine- 
tuning of timings is not required. It is suggested that the 
gap opening at which the satellites cease to form occurs 
during the period when satellites are trapped by the disk 
edge. This inferred gap opening time is consistent with 



that assumed in Section 13.21 In addition, we find that 
even if satellites are formed in the cold disk, the origin 
of the compositional gradient can be explained by colli- 
sional erosion of volatiles, as discussed in the companion 
paper. Furthermore, we can exclude the possibility that 
the Galilean satellites arc formed in the hot disk, that is, 
Jupiter did not open the gap in the protoplanetary disk 
when the circumplanctary disk was hot. 

4. DEPENDENCE ON PARAMETERS 

In the previous section, we presented -/V-body results of 
satellite accretion for our fiducial parameters and repro- 
duced the properties of the Galilean satellites (e.g., 2:1 
commensurability). There are, however, uncertainties in 
the parameters that we used. 

Gas surface de nsity (f g )'- Following 
iCanup fc Wardl (|2006D and iSasaki et all 
we assumed a = 5 x 10 3 and 
tq = 5 x 10 6 years, leading to a gas-starved 
disk with f g = 1 (Equation ©). This 
assumption of a small mass disk may be 
reasonable in the final stage of satellite 
formation because the inflow rate somehow 
decays, a nd gas surface d ensity decreases. 
Recently, iFuiii et all (|2011[ ) suggested a low 
magnetic Reynolds number and thus a low 
turbulent viscosity, which may suggest a 
more massive disk (/ g > 1 in Equation ([5])). 

Amount of solid materials (fd.m)'- We as- 
sumed /d,i n = 1 for the amount of solid 
materials, in which the gas-to-dust ratio in 
the inflowing gas is considered to be 100 ac- 
cording to solar metallicity. It is inferred 
that "metal" abundance in J upiter is high 
(e.g.. ISaumon fc Guillotll2004D : thus, the in- 
flowing gas would be more metal-rich (thus 
/dan > !)• This could be due to depletion of 
local gas by accretion onto Jupiter or that 
toward the Sun, leaving an excess of solid 
materials. However, dust grain abundance 
may also be reduced due to solids having 
been incorporated into planetesimals, which 
decreases /d,i n - Thus, depending on which 
effect is dominated, /d,i n can be larger or 
smaller than unity. 

Efficiency of type I migration (Cx)- In 
the fiducial model, the type I migration 
speed is set to the value derived from 
the linear calculation (Cx = 1). Several 
mechanisms for reducing migration efficiency 
have been discussed (e.g. JMasset et al1 F2006: 
iPaardekooper et ail 120101 ) . In fact, to repro- 
duce the observed distribu tions of exoplanets, 
CV < 0.1 (llda fc Lin|[200l . 

In this section, therefore, we first examine the mecha- 
nism through which the ./V-body results are altered by 
adopting different parameter values, and we attempt to 
link the model parameters to the physical properties of 
formed satellites, such as the resonant relationship, by 
using semianalytical arguments. This approach, in which 
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dependences of parameters on properties of satellites are 
investigated, is valid not only for setting limits on param- 
eters that reproduce the Galilean satellites but also for 
discussing satellite formation in general, including satel- 
lite formation around extrasolar giant planets. 

4.1. N-body Simulations with Various Parameters 

We present the results of iV-body simulations with var- 
ious parameters in Figures IU[a)-(e). The model parame- 
ters for each simulation are summarized in Table [2] The 
left panels of Figures [5] show orbital evolution of each 
representative run, while the right panels display aver- 
aged orbital distributions of the four innermost satellites 
at each time over three runs. An exception is Figure|6je) 
in which only the three innermost satellites are plotted. 
The times are selected after the time at which the inner- 
most satellite reaches the disk inner edge. 

Figure |U[a) shows the case of model 3, in which the 
gas surface density is 10 times higher than that adopted 
in model 1. Although the damping timescales of the 
eccentricity and the semimajor axis due to the gravita- 
tional drag become 10 times shorter, the trend agrees 
with that of model 1. Three to five satellites are lined up 
outside the disk inner edge, and they generally exhibit 
2:1 mean motion commensurabilitics. The critical mass 
for starting migration drops to ~ 2 x 10~ 5 Mp, and the 
masses of satellites trapped near the edge are comparable 
to those of the Galilean satellites (~ 5 x 10~ 5 Mp) before 
t = 6 x 10 6 Ir. After this time, the masses continue to 
increase as the satellites accrete inflowing solid materi- 
als. At t = 3 x 10 7 Ik, the mass of the satellites that 
are trapped by the edge become ~ 3 x 10~ 4 Mp, which 
is consistent with model 1 because the solid inflow rate 
is the same as that of model 1 (/dan = !)■ Despite the 
increase in the eccentricity damping force, the eccentric- 
ities of the satellites captured into resonances are nearly 
comparable to those in model 1, because the eccentrici- 
ties of the satellites exhibiti ng the eccentricity t rap are 
given as a function of t e /t a (Ogihara ct al. 2010h. 

Figures[6jb) and (c) show the results of models 4 and 5, 
respectively, in which the inflow rate of solid materials is 
changed. The critical mass for migration is proportional 

3/5 

to /, '. ; thus, the change of the inflow rate by a factor 
of 2 leads to a variation of the critical mass by a factor 
of 2 3 / 5 . The masses of satellites that are trapped by 
the edge at 3 x 10 7 T K are ~ 1.5 x 10 -4 M P (model 4) 
and ~ 6 x 10 -4 Mp (model 5). Although the masses 
of the migrating satellites change, altering the migration 
rate, the overall results are hardly affected. In particular, 
three to four satellites are captured into the 2:1 mean- 
motion resonance as is found in the fiducial model. 

FiguresEJd) and (e) show the results of models 6 and 7, 
in which the type I migration speed is reduced (model 6) 
and increased (model 7) from model 1 by a factor of 10. 
The critical mass is changes; however, the orbital config- 
uration in model 6, in which approximately four satellites 
in the 2:1 resonance are trapped by the edge, is approx- 
imately the same as in model 1. The masses of satellites 
trapped by the edge at 3 x 10 7 Ik are also comparable 
to those in model 1. However, the eccentricity trap is no 
longer effective in model 7, which results in orbital dis- 
tributions of satellites that differ significantly from those 
in other models. 




40 




Time (T ) Eccentricity 

Figure 6. Results of TV-body simulations of (a) model 3, (b) 
model 4, (c) model 5, (d) model 6, and (e) model 7. Left: Same as 
that described in FigureU with orbital evolution for representative 
runs. Right: Same as that described in Figure \3\ with averaged 
orbital distributions. 
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4.2. Resonances and Number of Satellites 

In this subsection, we discuss properties of remaining 
satellites in detail. To further discuss the resonant con- 
figuration and the total number of formed satellites over 
a wide parameter range, we calculate the orbital evo- 
lution of hypothetical three satellite systems. Here we 
show that capture into 2:1 resonance is robust for a rea- 
sonable range of parameters. The masses of satellites are 
set to 1CP 4 Mp. The gas surface density / g , which af- 
fects damping rates of e and a, and the migration speed 
C\ (hence t e /t a ) are treated as free parameters. The re- 
sults are summarized in Figure [7J Filled circles and filled 
squares in the figure indicate that all the three satellites 
exhibit the eccentricity trap and are captured into the 2:1 
and 3:2 resonances, respectively. Filled triangles repre- 
sent the cases in which the satellites exhibit the eccentric- 
ity trap but the inner pairs are in the 2:1 resonance and 
the outer pairs are in the 3:2 resonances. Inverted trian- 
gles represent inner pairs in the 3:2 resonance and outer 
pairs in the 2:1 resonances. Open squares represent tran- 
sitional cases from the filled triangles to the filled squares 
in which inner and outer pairs are temporarily captured 
into the 2:1 and 3:2 resonances, respectively, outside the 
disk edge; after a while, configurations of the inner pairs 
are converted to the 3:2 resonance when the innermost 
satellites move inside the edge. Crosses represent cases 
in which the eccentricity trap is not realized, although 
the satellites are captured into the 2:1 or 3:2 resonances. 
Asterisks show cases of satellites collisions. 
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Figure 7. Results of three body calculations for various values 
of / g and te/ta- Filled circles indicate that all three satellites are 
captured into the 2:1 resonance outside the disk inner edge. Solid, 
dashed, and dott ed li n es re prese nt se mianalytically estimated lim- 
its by Equations 1241 1 . d 3 ID - and 1321 1 . respectively. 



We find that the region (filled circles) in which satel- 
lites are captured into the 2:1 resonance outside the 
disk edge is vast. In fact, even with rapid migration 
(Ci ~ 2.5) and a high gas surface density (/ g ~ 10 2 ; 
£ g ~ 10 4 g cm -2 ), final configurations are identical. 
This region represents the necessary condition for the 
formation of the Galilean satellites, suggesting that a 
Galilean-like configuration in which the three satellites 
are in the 2:1 mean motion resonance can be formed with 
near certainty if the disk inner edge is sufficiently sharp. 
In some runs, we observe that the Laplace angle #5 li- 
brates about some fixed values. 

It is difficult to derive analytical formulae that rep- 
resent the filled circle region because resonant capture 



for equal-mass multiple systems has not been analyti- 
cally discussed. Therefore, we obtain semianalytical con- 
ditions by using simplified treatments. From Figure [7j 
we find that the following three conditions are necessary 
for capture into the 2:1 resonance outside the disk edge: 
(i) satellites should be trapped towing to the eccentric- 
ity trap to avoid falling onto the planet, (ii) eccentricity 
damping should not be excessively strong, and (iii) the 
migration speed should not be extremely rapid. Next we 
derive semianalytical formulae for each condition. 

Initially, the condition (i) is considered. For the ec- 
centricity trap to occur, a constraint on t e /t a (hence C\) 
is derived by an arg ument similar to that provided in 
lOgihara et al.l ()201ClO . If the edge torque, which is ex- 
erted on the innermost satellite because of the eccen- 
tricity damping force, balances migration torques of the 
outer satellites, then the satellites can be trapped by the 
edge. The trapping cond i tion i s obtained from Equa- 
tion (32) of lOgihara et~all ppToh : 



-A% el 1 
0.78 nT e ~ 2t~ 



1 



2t a 



>0, 



(22) 



where A c a {- 



-0.868 ) is the numerical coefficient 



(|Tanaka h Wardir2004f ). e is the eccentricity of the inner- 
most satellite, and n is the number of trapped satellites. 
The first term is the edge torque, the second term is the 
migration torque on the innermost satellite, and the third 
term is the resonant torque from outer satellites. Here, 
we assume that the satellites have nearly equal masses 
for simplicity. (Note that we used a more accurate ex- 
pr ession for the edge to rque than that in Equation (32) 
of lOgihara et~aT1 (I20T0I U 

From numerical calculations, eccentricities excited by 
resonances are well fitted by 



0.6 



1/2 



1/2 



(23) 



Substituting Equation (|23|) into Equation (|22j) . the fol- 
lowing constraint on t e /t a is obtained: 

<0.06Q)fAV. (24) 



n - 1 



ta 

where the second term of Equation (|22|) is ignored. The 
orbit of the innermost satellite partially enters the inner 
cavity, in which the satellite is not affected by the type I 
torque. In this case, the magnitude of type I migration 
torque on the innermost satellite should be smaller than 
l/2t a , and the third term dominates the second term 
(Ogih ara et al.ll2~010( ). The capture constraint by the ec- 
centricity trap is plotted with a solid line in Figure [7j 
This condition is consistent with the results of orbital 
calculations. 

Then, from condition (ii), we derive a constraint 
by comparing timescales of excitation and eccentricity 
damping. The eccentricity excitation by a single distant 
encounter of tw o satellites on approximately c ircular or- 
bits is given by (Hascgawa & Nakazawa 1990) 



A . I — ) 



(25) 



where b is the difference in the semimajor axes of the two 
satellites. Because the distant encounter occurs at every 
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synodic period 



2nr 



(26) 



syn " 360/2' 

the variation of eccentricity per unit time is expressed as 



de 
di 



For first-order resonances of j : j — 1, the orbital separa- 
tion is expressed as 



J 



3-1 



2/3 



1 



(28) 



then, Equation (|2~71) is reduced to 

M 



de 
dt 



1.7 x 10" 
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\20R P 



-3/2 



-3/2 



yr 



(29) 



where j = 2 is assumed in the final transformation. A 
necessary condition for capture into a resonance is that 
e-excitation due to distant encounters should be larger 
than e-damping due to the gas drag described in Equa- 
tion (JTUJl: 

de > -• (30) 



dt 



This leads to the following constraint on f g 
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(31) 
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From Figure [7J we find that the capture condition for 
the 2:1 resonance is well explained by / g < 700, in which 
e ~ 0.1 is used. This capture constraint is plotted with 
a dashed line in Figure [3 

Finally, from condition (iii), a necessary constraint 
condition on the migra tion speed is de rived. Accord- 
ing to both analytical ( Fricdland] 120011) and numerical 
(jlda et al.l [20001: lWvattll2003f ) studies of resonant trap- 
ping, capture probability depends on the mass and the 
semimajor axis. The dependence of the critical migration 
rate, which can be represented by Ci iCr it, for trapping is 
given as Ci iC1 it oc (M/Mp) -4 / 3 r2 -1 . Then, we adjust the 
coefficient by our numerical calculations. The numerical 
result shown in Figure [7] is fitted as 



Ci < 200/: 



M 



10" 4 Mp 



-4/3 



(32) 



which is the critical migration timescale of the three 
satellites for the 2:1 resonance, and it is plotted with 
a dotted line in Figure [JJ 

The enclosed shaded region in Figure [7] represents 
the condition for formation of Galilean-like satellites. 
The results of three satellite calculations and the semi- 
analytical arguments are consistent with the results of 
iV-body simulations. Figure [8] indicates the parameters 



used in our TV-body simulations on the gas surface den- 
sity (f g ) and the migration speed (Ci) plane. The shaded 
region corresponds to that in Figure [7j In case of mod- 
els 1-6, the parameters are in a good region for repro- 
ducing the number and the resonant relationships of the 
Galilean satellites, while the parameter values of model 7 
is far from the eccentricity trapping limit. 
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Figure 8. Model parameters adopted in iV-body simulations. The 
parameters of model 7 are outside the shaded region, which are 
derived by a semianalytical argument. 

We also compare our results with previous studies, in 
which the Hamiltqnian model is used in case of rapid 
migration ()Quillenll200l iMustill fc Wvattll2011[ ). These 
studies have considered the circular restricted three-body 
problem with a massive planet and a massless test par- 
ticle orbiting a central star to derive capture probability 
for first- and second-order resonances. Because the bod- 
ies with comparable masses are captured into resonances 
in our simulation and the effect of eccentricity damping is 
included, it is not possible to directly apply their predic- 
tions to our case. However, we find that their results are 
approximately consistent with our iV- body results and 
hence our semianalytical arguments. In IMustill fc Wvattl 
(|201 If) . capture probability is derived as a function of 

the generalized momentum J and the migration rate /3. 
Since J ~ 10 _1 and $ ~ 10~ 2 are derived in the case 
of model 1, it is determ ined from Figures 2 and 11 in 
IMustill fc Wvattl ([20111 ) that capture into the 2:1 reso- 
nance is highly likely. Therefore, our findings of a robust 
2:1 resonance formation are also confirmed by the H amil- 
tonian model. In addition, IMustill fc Wvattl ([201 1[) doc- 
umented the following trends, which is also observed in 
our calculations: capture probability decreases with a 
decrease in mass. If the migration rate increases by 10 3 
times of that in model 1, bodies tend to be captured into 
the closer 3:2 resonance. These features therefore agree 
with our semianalytical estimates. Further investigation 
of the resonant capture for comparable masses including 
eccentricity damping is expected in future work. 
Finally, we note that since the satellites are captured 
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into the 2:1 resonance, the orbital difference b is gener- 
al ly ~ 10— 15 fH, wh ich is larger than b = 5 rn assumed 
in fSasaki et al.l (|2010l ). In their calculations, small satel- 
lites are trapped between large satellites; as a result, or- 
bital differences between large satellites become ~ 15 rjj. 
Their assumption is valid for relatively fast migration 
cases (e.g., planet formation); therefore slight improve- 
ment is needed. This is also discussed in Section [5] 

4.3. Mass 

In this subsection, we discuss the satellite mass by 
using semianalytical estimates. As mentioned in Sec- 
tion 13.11 satellites grow according to the accretion 
timescale of Equation (p~9|) . Adding the exponential de- 
cay of the gas inflow, we can obtain the satellite mass as 
a function of time by integrating Equation (|19[) : 



M(t) 
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1 — exp 
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3/2 



(33) 



where Mq is the initial mass. The time dependence of 
r is neglected. Figure [§] plots the mass of the largest 
three satellites as a function of /d,in for 12 runs (mod- 
els 1, 3-5). The circles, triangles, and squares express the 
mass at 3 x 10 5 , 3 x 10 6 , and 3 x 10 7 Tk, respectively. 
The estimates of Equation (|3"3")) are also drawn with the 
dotted line (3 x 10 5 T K ), the dashed line (3 x 10 6 T K ), 
and the solid line (3 x 10 7 Ik), where r = 20 Rp is 
substituted. Because satellites migrate inward from the 
original locations, we neglect r-dependence. We find that 
the results are reasonably consistent with the estimates. 
However, it should be noted that although Ar = 10 rn is 
assumed in deriving Equation (|19j) . the orbital separation 
of satellites trapped by the edge can decrease to < 10 rn 
with an increase in mass, and Ar becomes independent 
of M. Thus, the power-law index in the right-hand-side 
of Equation (|3"3")l decreases to unity. In fact, the slope at 
3 x 10 7 Tk is clos e to th e linear dependence of /dan- 

ICanup fc Ward! (|2006l ) showed that the mass of satel- 
lites and the total mass in a system are both regulated 
by the balance between the supply of the inflowing ma- 
terials and orbital decay due to type I migration. The 
critical mass is derived in Equation (|2"0"|) . However, in 
our case in which the inner cavity exists, the total mass 
regulation does not work effectively unless the eccentric- 
ity trap is inhibited. For a case in which the number 
of trapped satellites by the edge is equal to the critical 
number for the eccentricity trap, the configuration is de- 
stroyed and some satellites are lost to the planet if an ad- 
ditional satellite migrates and is trapped in a resonance. 
In other cas es, the total mass c ontinues to increase. As 
discussed in lSasaki et al.l (|2010D . when the total mass of 
the trapped satellites exceeds the disk mass, the satel- 
lites may be released to the host planet, and mass regu- 
lation occurs. In addition , recent magnetohydr odynamic 
(MHD) simulations (e.g.. iRomanova et al.ll2008h suggest 
that the inner cavity may be buried by intermittent gas 
accretion onto the planet, and satellites trapped by the 
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Figure 9. Mass of the largest three satellites at 3 X 10 5 Tk (blue 
circles), 3 X 10 6 Tk (red triangles), and 3 X 10 7 Tk (black squares) 
for models 1 and 3 (/ g = 1), model 4 (/ g = 0.5), and model 5 
(/ g = 2). The dotted, dashed, and solid lines are estimates of 
Equation J33]l at 3 X 10 5 Tk, 3 X 10 6 T K , and 3 X 10 7 Tk, respec- 
tively. 

edge would episodically fall. These effects should be ex- 
amined in detail through MHD simulations. 

With an analytical estimate of the growth timescale, 
model parameters can be constrained for the formation of 
the Galilean satellites. To prevent these parameters from 
exceeding the masses of the Galilean satellites, M cr it < 
5 x 10 -5 Mp should be satisfied as a necessary condition, 
which leads to (Equation ([20])): 
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(34) 



where we neglect the impact erosion of ice. In addition, 
to be consistent with Callisto's partially differentiated 
state, an accretion time longer than ~ 5 x 10 5 years is 
implied to enable ?7ice/d,in < 1 at the time of Callisto's 
accretion (Equation (|T9|) V 



5. FORMATION OF SATURNIAN SYSTEM 



iSasaki et al.l (|2010D proposed that Jupiter opened a 
gap in the protoplanetary disk to halt its growth, and 
Saturn did not. The circum-Saturnian disk gradually 
dissipated on the global depletion timescale of the pro- 
toplanetary disk; therefore, in the final stage of satellite 
accretion around Saturn, the circum-Saturnian disk did 
not have an inner cavity because a low accretion rate 
would lead to a weak magnetic field, as discussed in Sec- 
tion [2l This condition is comparable with that adopted in 
ICanup fc Ward! (|2006f ): however, properties of the Satur- 
nian system are not commonly reproduced either by their 
Af-body simulations or our preliminary iV-body calcula- 
tions without the inner cavity. That is, the average num- 
ber of final satellites is ~ 3 — 4, while the number of large 
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satellites (> 10~ 5 Mp) in the current Saturnian system 
is one (Titan). Hence, the solid material distribution 
produced by iV-body calculations is also not consistent 
in the current Saturnian system. 

According to Equation p^|) . the growth timescale of 
satellites decreases with distance from the central planet; 
therefore, satellite seeds predominantly develop in the 
outer region (> 20 -Rp), and the satellites exhibit in- 
ward migration. If the migration rate significantly in- 
creases during migration, the migrating satellites rapidly 
fall onto the planet, and as a result, a single sa tellite al- 
ways exists in Titan's orbit. ISasaki et al.l (|2010f) assumed 
accelerated migration and showed that one or two satel- 
lites remain in the end; however, their method for gener- 
ating satellite seeds is not appropriate for satellite forma- 
tion. Because the migration timescale in protosatellite 
disks is proportional to r 1 / 2 (Equation (fTS"])). the migra- 
tion rate is not accelerated (da/dt — a/t a ~ a 1 / 2 ) during 
migration, which is ob served in TV-body simulations of 
iCanup fc Ward! (|2006l) and in this study. It should be 
noted that planets in a minimum-mass nebula, in which 
a steeper surface density profile (E g oc r 3 / 2 ) is consid- 
ered, undergo accelerated migration, because t a oc r 3 / 2 . 
Thus, we suggest in this section several possibilities to 
reconcile the inconsistency. 

In o ur ca lculations and also in those of lCanup fc Ward! 
(2006) and ISasaki et al.l (|2010f) . the radial dependence of 
the gas inflow is ignored, which results in a relatively 
gentle gradient of gas surface density (S g oc r -3 / 4 ). 
If the gas inflow flux is inversely proportional to r 
(Fin oc r _7 ;7 > 0), the gradient of gas surface den- 
sity becomes steep; re cent hydrodynamic simulations by 
iTanigawa et al.l (|2012l ) suggest that 7 ~ 1. In addition, 
photoev aporation can produce th e steep gradient of gas 
density ((Mitchell & Stewartll2"oH . The type I migration 
speed increases in the region close to the central planet 
so that the satellites that begin to undergo migration are 
rapidly lost to the planet. Our preliminary iV-body sim- 
ulations that include this effect suggest that the number 
of satellites decreases with an increase in 7. More de- 
tailed calculations are required, and it is also important 
to examine the radial and time dependence of the inflow 
flux by a high-resolution hydrodynamic simulation. 

Another possibility is that the solid density is locally 
increased around Titan's orbit (~ 20 Rp), leading to the 
preferential growth of satellites in the region. The radial 
distribution of solid inflow and satellitcsimal formation 
should be investigated. A detailed discussion of the Sat- 
urnian system formation will be presented in a separate 
paper. 

6. CONCLUSIONS 

We have investigated the formation of regular satel- 
lites around giant planets using direct iV-body simula- 
tions, including the effects of the eccentricity trap near 
the disk edge and the damping of the semimajor axis and 
eccentricity. Through these simulations, we confirmed 
the following scenarios of the final stage of a satellite for- 
mation in a steady accretion disk with an inner cavity 
and a gas/solid inflow from a circumstellar disk: 

1. Within the inflow region, protosatellites accrete 
from solid materials in their feeding zones. When 



the satellites grow to the critical mass for migra- 
tion, they begin to exhibit orbital decay. Because 
the width of the feeding zone increases with an in- 
crease in distance from the central planet and disk 
surface density has a relatively weak radial depen- 
dence (S g oc r~ 3 / 4 ), satellites in outer orbits grow 
larger than inner satellites. 

2. When the satellites reach the inner disk edge, they 
cease inward migration. The subsequently migrat- 
ing satellites are captured into mean motion reso- 
nances with the inner satellites. If the edge torque 
exerted on the innermost satellite balances the mi- 
gration torque, the satellites are trapped near the 
edge, which is called as the eccentricity trap. 

3. Satellites with the critical mass migrate inward to 
be trapped in successive resonances. If the number 
of trapped satellites exceeds the critical number for 
the eccentricity trap, the orbital configuration of 
the resonant convoy is disrupted, resulting in a loss 
of the inner satellites to the planet. That is, the 
total number of trapped satellites is self-regulated. 

4. The disk gas of circumplanetary protosatellite disks 
quickly (~ 10 3 years) vanishes because of trunca- 
tion of inflowing gas by the gap opening along the 
planet's orbit in the circumstellar disk, and reso- 
nant configuration is frozen at that time. Although 
the eccentricities of trapped satellites are excited 
by resonances, they are generally damped to < 
0.01 through long-term tidal dissipation within the 
satellites. 

We tracked the compositional evolution of satellites in 
iV-body simulations including the ice line beyond which 
water is condensed. As migration of icy satellites from 
the outer region ceases because of the eccentricity trap, 
inner satellites avoid contamination by water-rich mate- 
rials, resulting in a radial compositional gradient. This 
composition is described as an increase in icy components 
with distance. In addition, it is suggested that to repro- 
duce the compositional gradient, the gap opening needs 
to occur while the satellites are trapped by the edge via 
the eccentricity trap. An additional pathway for pro- 
ducing the compositional gradient is also suggested such 
that the volatile materials of inner satellites would be 
evaporated because of high-speed collisions during their 
accretion. 

Furthermore, we examine the relationship between the 
characteristics of the final satellites and the model pa- 
rameters by numerical simulations and semianalytical 
arguments to evaluate the probability of the Galilean- 
like satellite formation. The physical properties of the 
Galilean moons provide the following constraints on the 
model parameters. The number of satellites that can be 
trapped by the eccentricity trap restricts the ratio of the 
e-damping timescale versus the a-damping timescale: 

— < 0.05, (35) 

where n = 3 is substituted into Equation ((24)) . Equa- 
tion (|3"5"|) is given as (c s /vk) 2 < 0.01/Ci, which is usu- 
ally satisfied in the protosatellite disk with the temper- 
ature of Equation (0. The number of final satellites 
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also depends on the inflow region of the solid materi- 
als, although the actual inflow region of gas and dust 
for the circum- Jovi an orbit is pro bably consistent with 
our adopted theory (Machida 2009) . By considering cap- 
ture into resonances, constraints on the e-damping and 
a-damping timescales are derived. A necessary condition 
for capture into the 2:1 resonance is 



/ g < 700. 



(36) 



This parameter can be easily satisfied in the gas-starved 
disk (/ g = 1), and this condition implies that an actively- 
supplied disk with either a much lower viscosity and/or a 
more rapid inflow rate than considered in previous works 
can also produce 2:1 resonant orbits. An additional con- 
dition for satellites with masses of 5 x 10 -5 Mp 



lp is 



d < 500/g" 1 . 



(37) 



This parameter can also be satisfied because C\ < 1. We 
therefore find that formation of Galilean-like satellites is 
robust, because even with fast migration and a high gas 
surface density over three satellites, they are captured 
into the 2:1 mean motion resonance outside the disk inner 
edge. The parameters for growth and migration are also 
constrained by considering that the satellite mass does 
not exceed the critical mass for migration: 



Vicefd .1 



< 0.5. 



(38) 



It is also inferred that the solid inflow is decreased 
(?7ice/d.in < 1) at the time of Callisto's formation. In ad- 
dition, to prevent Callisto from being captured into mean 
motion resonances, the gap opening timing of Jupiter 
should be comparable to the timing of the completion 
of Callisto formation. This theory is consistent not only 
with the gap opening time, which was inferred from the 
discussion of composition, but also with the time as- 
sumed in our calculation. 

We finally provide theoretical predictions for character- 
istics of exomoons, which should be detected in the near 
future. If the disk inner edge width is sufficiently sharp, 
the eccentricity trap is likely to occur because the value 
of t e /t a is generally small even in the case that applies 
the full type I migration rate from the linear theory. The 
number of trapped satellites by the disk edge would be 
less than ~ 10 because the inflow region of solid material 
is probably con sistent with that in the adopted theory 
(Machida 2009). In the last stage of satellite accretion, 
the gas surface density would presumably be less than 
10 5 g cm -2 so that the damping timescales of eccentric- 
ity and semimajor axis are sufficiently long to enable the 
capture of satellites into the 2:1 mean-motion resonance, 



unless the mass of the satellite is significantly small. 

These properties of multiple satellite formation and 
capture into resonances imply that the satellites can be 
substantially affected by tidal heating. The satellites 
that have a sufficient mass to retain water in the HZ arc 
potentially habitable. Furthermore, according to tidal 
heating, their host planet is not required to be in the 
HZ; thus, the probability of forming Europa-like habit- 
able satellites, which have tidally heated oceans, would 
be high. From the perspective of tidally heated hab- 
itable moons, it is probable that satellites in relatively 
outer orbits retain water compared to close-in moons, 
which are strongly heated. Tidal heating is extremely 
important for habitable moons; therefore, further study 
of the evolution of satellites including tidal heating is nec- 
essary. By considering the migration efficiency and the 
accretion rate, the mass of normal satellite size would be 
~ 10 -4 Mp, which is similar to that of the Galilean satel- 
lites and Titan. It is possible for observable exomoons 
(> 0.2 M e ~6x 10~ 4 Mj) to form around giant planets 
more massive than Jupiter. Satellite systems with a large 
mass ratio of the satellite to the host planet are more 
stable for detecting exomoons (jKipping et al.|[2012f) . and 
when the factor f7i C e/d,in/Ci/g is larger than unity, larger 
satellites (> 10~ 4 M P ) are formed. 

We find that several satellites are formed locked in 
resonances under the assumpti on that the host p lanet 
opens a gap around its orbit. ISasaki et al.l ([2010D pro- 
posed that if the planet does not create a gap, satel- 
lites are unlikely to form in resonances, such as that ob- 
served in Saturnian satellites. Because a gap opening 
indicates that the planet may undergo type II migra- 
tion, it is expected that satellites systems in mean mo- 
tion resonances migrate inward to some extent. There- 
fore, planets that reside in closer orbits probably har- 
bor Jovian-syst em-like moons, although destabilization 
by tidal tor que (|Barnes fc O'BrieiJfeOOa or shrinkage of 
Hill sphere (|Namounil l2010') may also be important. In 
contrast, Saturnian-system-like moons orbit more distant 
planets. 
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APPENDIX 
ACTIVELY SUPPLIED DISK MODEL 



by 



On the basis of iCanup k. Wardl (|2002f ). the asymptotic disk gas surface density for a steady accretion disk is given 



0.55 



(Al) 



where is the outer edge of the diffused-out disk. The second and third terms in the bracket are correction terms 
due to outward diffusion. 

The disk is heated by luminosity from the central planet, viscous dissipation, and energy dissipation associated with 
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the difference between the free-fall energies of the incoming gas and a Keplerian orbit. We assume that heating is 
dominated by viscous dissipation; then, the photospheric temperature of the disk (Td) is determined by a balance 
between viscous heating and blackbody radiation from the photosphere: 

a SB T* ~ 9 -n^ g ~ °-^^tfF P (A2) 

where crge is the Stephan-Boltzman constant. When Fp = Mp/tq, the temperature is reduced to 

\MjJ \5xl0 6 jvj \20RpJ \Rj V ; 

The midplane temperature T is approximately given by T ~ (1 + St/S) 1 / 4 !^, where r is the optical depth. Because 
we are concerned with the disk temperature for evolution of the ice line, we assume that T ~ Td to avoid uncertainly 
in opacity. 

Thus, by adopting the alpha model for disk viscosity (v = ac s H), the gas surface density of the disk becomes 



5 x 10" 



5 x 10 6 yr 



-3/4 



(A5) 



GROWTH TIMES C ALE 

The growth rate of satellites is determined by the velocity dispersion a and the spatial density of satellitesimals p d . 
When velocity dispersion is smaller than the surface escape velocity of the satellite core, the accretion rate of the core 
with the mass M at r is 



M * nCR*p d „ * nCR 2 E d n K (= 2nC ( * 



Mp 



(Bl) 



where C is a numerical factor of 2-3 (jStewart fc Id a 2000) and Ed is the surface density of satellitesimals in the disk. 
Because R/r and Sd^ 2 /Mp are large compared to the corresponding values for planetary growth at ~ 1 AU around 
the Sun, the accretion rate is relatively high. Thus, the accretion timescale is 

M 
M 

c^y 1 ( s d 

2.5 



' acc.corc " 



(B2) 



;3 x 10 d 



50 g cm 



-2 



3 g cm 3 



1/3 



M 



10" 4 M F 



-1/3 



V M 3 



1/6 



(— 

\20Rp 



1/2 



Rp 



1/2 



e 

or 



(B3) 



where p and e are the internal density of the core and the RMS eccentricity of satellitesimals, respectively. Through 
numerical simulations, we find that Ed ~ 50 g cm -2 and e ~ 0.1. 

However, in the case in which the core accretion timescale T acCtCmc is smaller than the mass in flow timescale Ta c c. inflow , 
the growth timescale is controlled by the inflow flux of solid materials, as determined by iCanup fc Wardl ([2006). 
Satellites accrete materials across an annulus of width Ar and grow in the timescale as 



_M 

^acc, inflow — ~ 

M 



fM 



F m 2nrAr 

l5„-l f-l 



20i? P 



M 



10" 4 Mp 



2/3 



yr, 



(B4) 
(B5) 



where /(= 100/r?i ce .f d.in) is the gas-to-solid mass ratio in the inflow, and Ar = 10 rn is used. This derivation differs 
slightly from that of ICanup fc W ard (20Q6|), in which Ar ~ 2er is assumed. However, the timescale is approximately 
consistent. Obviously, r acc ^ n fl ow is longer than r acCjCorc ; thus, r acCi i n fl ow can be considered as the growth timescale of 
the satellite. These estimates are in good agreement with numerical results shown in the text. 



The angels are defined as 



RESONANT ANGLE 



h = Ai — 2A2 + wi, 



(CI) 
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@5 — 62 — $3 — Ai — 3A2 + 2A3, 



62 = Ai — 2A2 + 

#3 = A2 — 2A 3 + zu 2 , 

^4 = A2 — 2 A3 + TU3, 



(C2) 
(C3) 
(C4) 
(C5) 



where \ and Wi are mean orbital longitudes and longitudes of the pericenter, respectively. The subscripts i = 1, 2, and 
3 refer to the innermost satellite (body 1), the second satellite (body 2), and the third satellite (body 3), respectively. 
These angles represent displacements of longitude of conjunction from pericenters; thus, 6\ = 0°, 9 2 — 180° indicates 
that conjunctions of body 1 and body 2 occur when body 1 is near the apocenter and body 2 is near the pericenter. 
If the angles librate about some fixed values, it is usually considered that they are in mean motion resonances. In 
addition, the Laplace resonance is characterized by the libration of 65. In this case, conjunctions of body 1 and body 2 
drift at the same rate as those of body 2 and body 3 so that triple conjunction is impossible. 

For the Galilean satellites, Q x an d #2 librate about 0° and 180°, respectively, indicating that the conjunction longitude 
between Io and Europa is locked to Io's pericenter and Europa's apocenter. 6 3 librates about 0°, while 6*4 circulates 
through 360°, which indicates that the Europa-Ganymede conjunction is locked to Europa's pericenter but to neither 
apse of Ganymede. Thus, #5 obviously librates about 180°, indicating a Laplace relationship. 
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